home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / Utilities / Quick3 / MonotonicFit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-23  |  4.9 KB  |  154 lines  |  [TEXT/KAHL]
  1. /*
  2. MonotonicFit.c
  3. Copyright (c) 1979,1989,1990 Denis G. Pelli 
  4.  
  5. FIND MAXIMUM LIKELIHOOD MONOTONIC FIT (MLMF)
  6. Given the number of right and wrong responses at each of several
  7. contrasts, this routine finds the maximum likelihood monotonic
  8. psychometric function that could have given rise to the data.
  9.  
  10. The original dataRecord is OVERWRITTEN and replaced by the "collapsed"
  11. data which represent the maximum likelihood monotonic fit.
  12.  
  13. The log likelihood and degrees of freedom are returned, using the
  14. user-supplied points. The degrees of freedom is equal to the number of
  15. points in the maximum likelihood monotonic fit (after collapsing), 
  16. minus 1 if both limiting probabilities (0 and 1) are included.
  17.  
  18. THEOREM: The MLMF when two successive points have a nonmonotonic
  19. proportion correct is equivalent to the MLMF to collapsed data where
  20. the two points have been combined into one. "Equivalent" means that
  21. the two MLMFs will be equal at all points other than the two in
  22. question, and the original MLMF will have the same value at both of
  23. the two nonmonotonic points as the equiv. MLMF will have at the
  24. combined point. Repeated application of the theorem eventually gives
  25. collapsed data with monotonically increasing proportion correct. The
  26. MLMF is equal to these proportion correct.
  27.  
  28. The principal use of the MLMF routine is to deal with the problem of assessing
  29. the goodness of a function fit (e.g. Weibull) to data
  30. collected by a staircase that has put different numbers of trials at
  31. each contrast. Clearly the unconstrained fit is sort of silly for this case, 
  32. especially when there are many
  33. contrasts with only single trials, and it's not clear how to count
  34. its degrees of freedom. (One contrast may have 200 trials, another may
  35. have one trial. Do they each count as a degree of freedom? No good answer.)
  36. The MLMF finds the maximum likelihood monotonic
  37. fit, (and along the way determines the number of degrees of freedom), 
  38. and in practice makes a good null hypothesis for
  39. assessing the Weibull fit. What's cute about MLMF is that it is fast.
  40. In time proportional to the number of contrasts it computes the 
  41. monotonic fit directly. I worked this out for myself while in graduate 
  42. school, but Misha Pavel tells me that there is actually a whole literature
  43. on this topic.
  44. HISTORY:
  45. 29-JUL-79    dgp    wrote it in FORTRAN
  46. 8/9/89         dgp    translated it to C, but but didn't test it.
  47. 4/5/90        dgp    debugged.
  48. 4/7/90        dgp    added qsort().
  49. 10/29/90    dgp    tidied up the comments.
  50. */
  51. #include "VideoToolbox.h"
  52. #include "Quick3.h"
  53.  
  54. void MonotonicFit(dataRecord *dataPtr,double *logLikelihoodPtr,int *degreesOfFreedomPtr)
  55. {
  56.     enum{END=-1,DELETED=-2};
  57.     int head;
  58.     long y[MAX_CONTRASTS];        /* long to allow huge number of trials per point */
  59.     long t[MAX_CONTRASTS];
  60.     long n;
  61.     double p[MAX_CONTRASTS];    /* the values of the MLMF psychometric function. */
  62.     int next[MAX_CONTRASTS];
  63.     int flag;
  64.     int i,j,jj;
  65.  
  66.     /*
  67.     Sort the contrastRecords in order of increasing contrast
  68.     & merge data at equal contrasts.
  69.     */
  70.     SortAndMergeContrasts(dataPtr);
  71.  
  72.     /* Load up the arrays. next[] is a linked list. */
  73.     head=0;
  74.     for(i=0;i<dataPtr->contrasts;i++){
  75.         y[i]=dataPtr->c[i].correct;
  76.         t[i]=dataPtr->c[i].trials;
  77.         p[i]=y[i]/(double)t[i];
  78.         next[i]=i+1;
  79.     }
  80.     next[i-1]=END;
  81.     /*
  82.     Assume the data are in order of increasing contrast.
  83.     Now repetitively apply the theorem.
  84.     Each time a nonmonotonic pair of points is found,
  85.     they are replaced by a single point with the sum of the
  86.     data at the two points.
  87.     Successive passes are made through the linked list until
  88.     no nonmonotonicities are found (i.e. flag==FALSE).
  89.     */
  90.     flag=TRUE;
  91.     while(flag){
  92.         flag=FALSE;
  93.         j=head;
  94.         for(i=0;i<dataPtr->contrasts;i++){
  95.             jj=next[j];
  96.             if(jj==END)break;
  97.             if(t[j]==0 || t[jj]==0 || p[j] > p[jj]) {
  98.                 /* these two points are nonmonotonic, so merge them */
  99.                 /* Points with zero trials should be merged as well */
  100.                 flag=TRUE;                            
  101.                 y[j] += y[jj];
  102.                 t[j] += t[jj];
  103.                 p[j]=y[j]/(double)t[j];
  104.                 next[j]=next[jj];
  105.                 next[jj]=DELETED;                /* mark the deleted point */
  106.                 jj=j;
  107.             }
  108.             j=jj;
  109.         }
  110.     }
  111.  
  112.     /* 
  113.     Assign probabilities to the deleted points.
  114.     */
  115.     for(i=1;i<dataPtr->contrasts;i++){
  116.         if(next[i]==DELETED) p[i]=p[i-1];
  117.     }
  118.  
  119.     /* 
  120.     Count the number of degrees of freedom.
  121.     */
  122.     *degreesOfFreedomPtr=0;
  123.     for(i=0;i<dataPtr->contrasts;i++){
  124.         if(next[i]!=DELETED) (*degreesOfFreedomPtr)++;
  125.     }
  126.     if(p[0]==0.0 && p[dataPtr->contrasts-1] == 1.0) (*degreesOfFreedomPtr)--;
  127.  
  128.     /*
  129.     Compute the log likelihood
  130.     */
  131.     *logLikelihoodPtr=0.0;
  132.     for(i=0;i<dataPtr->contrasts;i++){
  133.         n=dataPtr->c[i].correct;
  134.         if(n>0) *logLikelihoodPtr +=n*log(p[i]);
  135.         n=dataPtr->c[i].trials-dataPtr->c[i].correct;
  136.         if(n>0) *logLikelihoodPtr +=n*log(1.0-p[i]);
  137.     }
  138.  
  139.     /*
  140.     To return the results, OVERWRITE the dataRecord with
  141.     the maximum likelihood monotonic fit, where p=correct/trials.
  142.     */
  143.     for(i=0;i<dataPtr->contrasts;i++){
  144.         if(next[i]!=DELETED) {
  145.             dataPtr->c[i].correct=y[i];
  146.             dataPtr->c[i].trials=t[i];
  147.         }else{
  148.             dataPtr->c[i].correct=dataPtr->c[i-1].correct;
  149.             dataPtr->c[i].trials=dataPtr->c[i-1].trials;
  150.         }
  151.     }
  152. }
  153.  
  154.